---
title: "7-7 mMRC breathlessness scale"
description: |
Analyses of the odified Medical Research Council breathlessness scale at 28 days.
author:
- name: James Totterdell
affiliation: University of Sydney
- name: Rob Mahar
affiliation: University of Melbourne
date: today
toc-depth: 5
editor_options:
chunk_output_type: console
---
```{r}
#| label: pkgs
#| code-summary: Load packages
library (ASCOTr)
library (tidyverse)
library (lubridate)
library (kableExtra)
library (patchwork)
library (cmdstanr)
library (posterior)
library (bayesplot)
library (ggdist)
library (lme4)
library (broom)
library (broom.mixed)
library (bayestestR)
theme_set (theme_classic (base_size = 10 , base_family = "Palatino" ) +
theme (panel.grid = element_blank (),
strip.background = element_blank ()))
bayesplot_theme_set (theme_set (theme_classic (base_size = 10 , base_family = "Palatino" ) +
theme (panel.grid = element_blank (),
strip.background = element_blank ())))
color_scheme_set ("red" )
options (digits = 4 )
```
```{r}
#| label: analysis-sets
#| code-summary: Prepare datasets
all_data <- readRDS (file.path (ASCOT_DATA, "all_data.rds" ))
breath_labels <- c ("Only breathless with strenuous exercise" ,
"Short of breath when hurrying on level ground or walking up a slight hill" ,
"Walks slower than most people of the same age because of breathlessness on level" ,
"Stops for breath after walking about 100 metres or after a few minutes on level ground" ,
"Too breathless to leave the house, or breathless when dressing or undressing" )
breath_labels_short <- c ("With exsercise" ,
"Up a slight hill" ,
"Slow for age" ,
"After 100 metres" ,
"Can't leave house" )
# FAS-ITT
fas_itt_dat <- ASCOTr::: make_fas_itt_set (all_data) %>%
mutate (out_mmrc_lab = fct_explicit_na (
factor (out_mmrc_scale,
levels = c (1 : 5 , 0 ),
labels = c (breath_labels_short, "Not asked" ))))
fas_itt_nona_dat <- fas_itt_dat %>%
filter (! is.na (out_mmrc_scale))
# ACS-ITT
acs_itt_dat <- ASCOTr::: make_acs_itt_set (all_data) %>%
mutate (out_mmrc_lab = fct_explicit_na (
factor (out_mmrc_scale,
levels = c (1 : 5 , 0 ),
labels = c (breath_labels_short, "Not asked" ))))
acs_itt_nona_dat <- acs_itt_dat %>%
filter (! is.na (out_mmrc_scale))
# AVS-ITT
avs_itt_dat <- ASCOTr::: make_avs_itt_set (all_data) %>%
mutate (out_mmrc_lab = fct_explicit_na (
factor (out_mmrc_scale,
levels = c (1 : 5 , 0 ),
labels = c (breath_labels_short, "Not asked" ))))
avs_itt_nona_dat <- avs_itt_dat %>%
filter (! is.na (out_mmrc_scale))
```
```{r}
#| label: models
#| code-summary: Load models
ordmod0 <- compile_cmdstanr_mod (
file.path ("ordinal" , "logistic_cumulative" ), dir = "stan" )
ordmod <- compile_cmdstanr_mod (
file.path ("ordinal" , "logistic_cumulative_epoch_site" ), dir = "stan" )
logistic <- compile_cmdstanr_mod (
file.path ("binary" , "logistic_site_epoch" ), dir = "stan" )
```
```{r}
#| label: helpers
#| code-summary: Helper functions
make_summary_table_anticoagulation <- function (dat, format = "html" ) {
tdat <- dat %>%
group_by (
CAssignment = factor (CAssignment, levels = c ("C0" , "C1" , "C2" , "C3" , "C4" ), labels = intervention_labels2 ()$ CAssignment)) %>%
summarise (
Patients = n (),
Known = sum (! is.na (out_mmrc_scale)),
` With exercise ` = sprintf ("%i (%.0f%%)" ,
sum (out_mmrc_scale == 1 , na.rm = TRUE ),
100 * mean (out_mmrc_scale == 1 , na.rm = TRUE )),
` Up a slight hill ` = sprintf ("%i (%.0f%%)" ,
sum (out_mmrc_scale == 2 , na.rm = TRUE ),
100 * mean (out_mmrc_scale == 2 , na.rm = TRUE )),
` Slow for age ` = sprintf ("%i (%.0f%%)" ,
sum (out_mmrc_scale == 3 , na.rm = TRUE ),
100 * mean (out_mmrc_scale == 3 , na.rm = TRUE )),
` After 100 metres ` = sprintf ("%i (%.0f%%)" ,
sum (out_mmrc_scale == 4 , na.rm = TRUE ),
100 * mean (out_mmrc_scale == 4 , na.rm = TRUE )),
` Can't leave house ` = sprintf ("%i (%.0f%%)" ,
sum (out_mmrc_scale == 5 , na.rm = TRUE ),
100 * mean (out_mmrc_scale == 5 , na.rm = TRUE )),
` mMRC scale, Median (Q1, Q3) ` = sprintf (
"%.0f (%.0f, %.0f)" ,
median (out_mmrc_scale, na.rm = T),
quantile (out_mmrc_scale, 0.25 , na.rm = TRUE ),
quantile (out_mmrc_scale, 0.75 , na.rm = TRUE ))
) %>%
bind_rows (
dat %>%
group_by (CAssignment = "Overall" ) %>%
summarise (
Patients = n (),
Known = sum (! is.na (out_mmrc_scale)),
` With exercise ` = sprintf ("%i (%.0f%%)" ,
sum (out_mmrc_scale == 1 , na.rm = TRUE ),
100 * mean (out_mmrc_scale == 1 , na.rm = TRUE )),
` Up a slight hill ` = sprintf ("%i (%.0f%%)" ,
sum (out_mmrc_scale == 2 , na.rm = TRUE ),
100 * mean (out_mmrc_scale == 2 , na.rm = TRUE )),
` Slow for age ` = sprintf ("%i (%.0f%%)" ,
sum (out_mmrc_scale == 3 , na.rm = TRUE ),
100 * mean (out_mmrc_scale == 3 , na.rm = TRUE )),
` After 100 metres ` = sprintf ("%i (%.0f%%)" ,
sum (out_mmrc_scale == 4 , na.rm = TRUE ),
100 * mean (out_mmrc_scale == 4 , na.rm = TRUE )),
` Can't leave house ` = sprintf ("%i (%.0f%%)" ,
sum (out_mmrc_scale == 5 , na.rm = TRUE ),
100 * mean (out_mmrc_scale == 5 , na.rm = TRUE )),
` mMRC scale, Median (Q1, Q3) ` = sprintf (
"%.0f (%.0f, %.0f)" ,
median (out_mmrc_scale, na.rm = T),
quantile (out_mmrc_scale, 0.25 , na.rm = TRUE ),
quantile (out_mmrc_scale, 0.75 , na.rm = TRUE ))
)
) %>%
rename (` Anticoagulation \n intervention ` = CAssignment)
kable (
tdat,
format = format,
align = "lrrrrr" ,
booktabs = TRUE ,
linesep = ""
) %>%
kable_styling (
font_size = 9 ,
latex_options = "HOLD_position"
) %>%
row_spec (nrow (tdat), bold = T)
}
make_summary_table_antiviral <- function (dat, format = "html" ) {
tdat <- dat %>%
group_by (
AAssignment = factor (AAssignment, levels = c ("A0" , "A1" , "A2" ), labels = intervention_labels2 ()$ AAssignment)) %>%
summarise (
Patients = n (),
Known = sum (! is.na (out_mmrc_scale)),
` With exercise ` = sprintf ("%i (%.0f%%)" ,
sum (out_mmrc_scale == 1 , na.rm = TRUE ),
100 * mean (out_mmrc_scale == 1 , na.rm = TRUE )),
` Up a slight hill ` = sprintf ("%i (%.0f%%)" ,
sum (out_mmrc_scale == 2 , na.rm = TRUE ),
100 * mean (out_mmrc_scale == 2 , na.rm = TRUE )),
` Slow for age ` = sprintf ("%i (%.0f%%)" ,
sum (out_mmrc_scale == 3 , na.rm = TRUE ),
100 * mean (out_mmrc_scale == 3 , na.rm = TRUE )),
` After 100 metres ` = sprintf ("%i (%.0f%%)" ,
sum (out_mmrc_scale == 4 , na.rm = TRUE ),
100 * mean (out_mmrc_scale == 4 , na.rm = TRUE )),
` Can't leave house ` = sprintf ("%i (%.0f%%)" ,
sum (out_mmrc_scale == 5 , na.rm = TRUE ),
100 * mean (out_mmrc_scale == 5 , na.rm = TRUE )),
` mMRC scale, Median (Q1, Q3) ` = sprintf (
"%.0f (%.0f, %.0f)" ,
median (out_mmrc_scale, na.rm = T),
quantile (out_mmrc_scale, 0.25 , na.rm = TRUE ),
quantile (out_mmrc_scale, 0.75 , na.rm = TRUE ))
) %>%
bind_rows (
dat %>%
group_by (AAssignment = "Overall" ) %>%
summarise (
Patients = n (),
Known = sum (! is.na (out_mmrc_scale)),
` With exercise ` = sprintf ("%i (%.0f%%)" ,
sum (out_mmrc_scale == 1 , na.rm = TRUE ),
100 * mean (out_mmrc_scale == 1 , na.rm = TRUE )),
` Up a slight hill ` = sprintf ("%i (%.0f%%)" ,
sum (out_mmrc_scale == 2 , na.rm = TRUE ),
100 * mean (out_mmrc_scale == 2 , na.rm = TRUE )),
` Slow for age ` = sprintf ("%i (%.0f%%)" ,
sum (out_mmrc_scale == 3 , na.rm = TRUE ),
100 * mean (out_mmrc_scale == 3 , na.rm = TRUE )),
` After 100 metres ` = sprintf ("%i (%.0f%%)" ,
sum (out_mmrc_scale == 4 , na.rm = TRUE ),
100 * mean (out_mmrc_scale == 4 , na.rm = TRUE )),
` Can't leave house ` = sprintf ("%i (%.0f%%)" ,
sum (out_mmrc_scale == 5 , na.rm = TRUE ),
100 * mean (out_mmrc_scale == 5 , na.rm = TRUE )),
` mMRC scale, Median (Q1, Q3) ` = sprintf (
"%.0f (%.0f, %.0f)" ,
median (out_mmrc_scale, na.rm = T),
quantile (out_mmrc_scale, 0.25 , na.rm = TRUE ),
quantile (out_mmrc_scale, 0.75 , na.rm = TRUE ))
)
) %>%
rename (` Antiviral \n intervention ` = AAssignment)
kable (
tdat,
format = format,
align = "lrrrrr" ,
booktabs = TRUE ,
linesep = ""
) %>%
kable_styling (
font_size = 9 ,
latex_options = "HOLD_position"
) %>%
row_spec (nrow (tdat), bold = T)
}
```
# Descriptive {#descriptive}
## Anticoagulation
```{r}
#| label: tbl-7-7-summary-anticoagulation
#| code-summary: Summary of mMRC outcome by arm
#| tbl-cap: Summary of mMRC scale at day 28 by treatment group, anticoagulation domain, FAS-ITT.
save_tex_table (
make_summary_table_anticoagulation (acs_itt_dat, "latex" ),
file.path ("outcomes" , "secondary" , "7-7-anticoagulation-summary" ))
make_summary_table_anticoagulation (acs_itt_dat)
```
```{r}
p1 <- fas_itt_nona_dat %>%
group_by (CAssignment = factor (CAssignment, labels = str_replace (intervention_labels ()$ CAssignment, "<br>" , " \n " ))) %>%
summarise (n = n ()) %>%
ggplot (., aes (CAssignment, n)) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Anticoagulation intervention" )
p2 <- fas_itt_nona_dat %>%
filter (! is.na (out_mmrc_scale)) %>%
dplyr:: count (
CAssignment = factor (
CAssignment,
labels = str_replace (intervention_labels ()$ CAssignment, "<br>" , " \n " )),
out_mmrc_scale = ordered (as.integer (out_mmrc_scale), levels = 0 : 5 )
) %>%
group_by (CAssignment) %>%
mutate (p = n / sum (n)) %>%
ggplot (., aes (CAssignment, p)) +
geom_col (aes (fill = out_mmrc_scale)) +
scale_fill_viridis_d (option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = F)) +
labs (fill = "mMRCbs" , y = "Cumulative proportion" , x = "Anticoagulation intervention" )
p <- p1 | p2
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "7-7-anticoagulation.pdf" ), p2, height = 3 , width = 6 )
p2
```
## Antiviral
```{r}
#| label: tbl-7-7-summary-antiviral
#| code-summary: Summary of mMRC outcome by arm
#| tbl-cap: Summary of mMRC scale at day 28 by treatment group, antiviral domain, FAS-ITT.
save_tex_table (
make_summary_table_antiviral (avs_itt_dat, "latex" ),
file.path ("outcomes" , "secondary" , "7-7-antiviral-summary" ))
make_summary_table_antiviral (avs_itt_dat)
```
```{r}
p1 <- fas_itt_nona_dat %>%
group_by (AAssignment = factor (AAssignment, labels = str_replace (intervention_labels ()$ AAssignment, "<br>" , " \n " ))) %>%
summarise (n = n ()) %>%
ggplot (., aes (AAssignment, n)) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Antiviral intervention" )
p2 <- fas_itt_nona_dat %>%
filter (! is.na (out_mmrc_scale)) %>%
dplyr:: count (
AAssignment = factor (
AAssignment,
labels = str_replace (intervention_labels ()$ AAssignment, "<br>" , " \n " )),
out_mmrc_scale = ordered (as.integer (out_mmrc_scale), levels = 0 : 5 )
) %>%
group_by (AAssignment) %>%
mutate (p = n / sum (n)) %>%
ggplot (., aes (AAssignment, p)) +
geom_col (aes (fill = out_mmrc_scale)) +
scale_fill_viridis_d (option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = F)) +
labs (fill = "mMRCbs" , y = "Cumulative proportion" , x = "Antiviral intervention" )
p <- p1 | p2
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "7-7-antiviral.pdf" ), p2, height = 3 , width = 6 )
p
```
## Age
```{r}
#| label: fig-7-7-age
#| fig-cap: |
#| Distribution of mMRC scale at day 28 by age group, FAS-ITT
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
dplyr:: count (
agegrp = cut (AgeAtEntry, c (18 , seq (25 , 75 , 5 ), 100 ), include.lowest = T, right = F),
mrc = ordered (out_mmrc_scale, levels = 0 : 5 ),
.drop = F) %>%
group_by (agegrp) %>%
mutate (p = n / sum (n))
pdat2 <- pdat %>%
group_by (agegrp) %>%
summarise (n = sum (n))
p1 <- ggplot (pdat2, aes (agegrp, n)) +
geom_col (colour = "grey40" , fill = "grey40" ) +
geom_vline (xintercept = 60 , linetype = 2 ) +
labs (y = "Number of \n participants" ,
x = "Age at entry" ) +
geom_vline (xintercept = 8.5 , linetype = 2 ) +
theme (axis.text.x = element_text (hjust = 1 , angle = 45 ))
p2 <- ggplot (pdat, aes (agegrp, p)) +
geom_col (aes (fill = mrc)) +
labs (x = "Age" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("mMRCbs" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = F)) +
geom_vline (xintercept = 8.5 , linetype = 2 ) +
theme (axis.text.x = element_text (hjust = 1 , angle = 45 )) +
theme (legend.key.size = unit (0.5 , "lines" ))
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-7-age.pdf" ), p, height = 2.5 , width = 6 )
p
```
## Sex
```{r}
#| label: fig-7-7-sex
#| fig-cap: |
#| Distribution of mMRC scale at day 28 by sex, FAS-ITT
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
dplyr:: count (
Sex,
mrc = ordered (out_mmrc_scale, levels = 0 : 5 ),
.drop = F) %>%
group_by (Sex) %>%
mutate (p = n / sum (n))
pdat2 <- pdat %>%
group_by (Sex) %>%
summarise (n = sum (n))
p1 <- ggplot (pdat2, aes (Sex, n)) +
geom_col (colour = "grey40" , fill = "grey40" ) +
labs (y = "Number of \n participants" ,
x = "Sex" )
p2 <- ggplot (pdat, aes (Sex, p)) +
geom_col (aes (fill = mrc)) +
labs (x = "Sex" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("mMRCbs" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = F)) +
theme (legend.key.size = unit (0.5 , "lines" ))
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-7-sex.pdf" ), p, height = 2.5 , width = 6 )
p
```
## Oxygen
```{r}
#| label: fig-7-7-oxygen
#| fig-cap: |
#| Distribution of mMRC scale at day 28 by supplemental oxygen requirement, FAS-ITT
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
dplyr:: count (
supp_oxy2 = fct_explicit_na (factor (supp_oxy2, levels = 0 : 1 , labels = c ("Not required" , "Required" ))),
mrc = ordered (out_mmrc_scale, levels = 0 : 5 ),
.drop = F) %>%
group_by (supp_oxy2) %>%
mutate (p = n / sum (n))
pdat2 <- pdat %>%
group_by (supp_oxy2) %>%
summarise (n = sum (n))
p1 <- ggplot (pdat2, aes (supp_oxy2, n)) +
geom_col (colour = "grey40" , fill = "grey40" ) +
labs (y = "Number of \n participants" ,
x = "Supplemental oxygen" )
p2 <- ggplot (pdat, aes (supp_oxy2, p)) +
geom_col (aes (fill = mrc)) +
labs (x = "Supplemental oxygen" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("mMRCbs" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = F)) +
theme (legend.key.size = unit (0.5 , "lines" ))
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-7-oxygen.pdf" ), p, height = 2.5 , width = 6 )
p
```
## Country
```{r}
#| label: fig-7-7-country
#| fig-cap: Distribution of mMRC scale at day 28 by country, FAS-ITT.
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
dplyr:: count (
Country = factor (PT_CountryName, levels = c ("India" , "Australia" , "Nepal" , "New Zealand" ),
labels = c ("India" , "Australia" , "Nepal" , "New \n Zealand" )),
mrc = ordered (out_mmrc_scale, levels = 0 : 5 ),
.drop = F) %>%
group_by (Country) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
pdat2 <- pdat %>%
group_by (Country) %>%
summarise (n = sum (n))
p1 <- ggplot (pdat2, aes (Country, n)) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Country of enrolment" )
p2 <- ggplot (pdat, aes (Country, p)) +
geom_col (aes (fill = mrc)) +
labs (x = "Country" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("mMRCbs" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = F)) +
theme (legend.key.size = unit (0.5 , "lines" ))
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-7-country.pdf" ), p, height = 2.5 , width = 6 )
p
```
## Site
```{r}
#| label: fig-7-7-site
#| fig-cap: |
#| Distribution of mMRC scale by study site, FAS-ITT.
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
dplyr:: count (
Country = factor (PT_CountryName, levels = c ("India" , "Australia" , "Nepal" , "New Zealand" ),
labels = c ("India" , "Australia" , "Nepal" , "New \n Zealand" )),
Site = fct_infreq (Location),
mrc = ordered (out_mmrc_scale, levels = 0 : 5 )) %>%
complete (mrc = ordered (0 : 5 ), nesting (Country, Site), fill = list (n = 0 )) %>%
group_by (Country, Site) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup () %>%
mutate (
Country = droplevels (Country),
Site = droplevels (Site)
)
pdat2 <- pdat %>%
group_by (Country, Site) %>%
summarise (n = sum (n)) %>%
ungroup ()
p1 <- ggplot (pdat2, aes (Site, n)) +
facet_grid ( ~ Country, scales = "free_x" , space = "free_x" ) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "" ) +
theme (axis.text.x = element_text (angle = 45 , hjust = 1 ),
panel.border = element_rect (fill = NA ))
p2 <- ggplot (pdat, aes (Site, p)) +
facet_grid ( ~ Country, scales = "free_x" , space = "free_x" ) +
geom_col (aes (fill = mrc)) +
labs (x = "Site" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("mMRCbs" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = F, ncol = 1 )) +
theme (axis.text.x = element_text (hjust = 1 , angle = 45 )) +
theme (legend.key.size = unit (0.25 , "lines" ))
p <- p1 / p2 +
plot_layout (guides = 'collect' )
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-7-country-site.pdf" ), p, height = 4 , width = 6.25 )
p
```
## Calendar Time
```{r}
#| label: fig-7-7-calendar
#| fig-cap: |
#| Distribution of mMRC scale by calendar time, FAS-ITT.
#| fig-height: 4
pdat <- fas_itt_nona_dat %>%
filter (! is.na (out_mmrc_scale)) %>%
dplyr:: count (
yr = year (RandDate), mth = month (RandDate),
mrc = ordered (out_mmrc_scale, levels = 0 : 5 ),
.drop = F) %>%
group_by (yr, mth) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p1 <- pdat %>%
group_by (yr, mth) %>%
summarise (n = sum (n)) %>%
ggplot (., aes (mth, n)) +
facet_grid ( ~ yr, drop = T, scales = "free_x" , space = "free" ) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Calendar date (month of year)" ) +
scale_x_continuous (breaks = 1 : 12 )
p2 <- ggplot (pdat, aes (mth, p)) +
facet_grid ( ~ yr, drop = T, scales = "free_x" , space = "free" ) +
geom_col (aes (fill = mrc)) +
labs (x = "Calendar date (month of year)" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("mMRCbs" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = F)) +
theme (legend.key.size = unit (0.5 , "lines" )) +
scale_x_continuous (breaks = 1 : 12 )
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-7-calendar-time.pdf" ), p, height = 2 , width = 6 )
p
```
## Sample Cumulative Logits
Proportional odds looks reasonable for all logits.
```{r}
#| code-summary: Plot sample cumulative logits
#| fig-cap: Inspect sample cumulative logits, anticoagulation.
trt_counts <- fas_itt_nona_dat %>%
dplyr:: count (CAssignment, out_mmrc_scale) %>%
complete (CAssignment, out_mmrc_scale, fill = list (n = 0 )) %>%
group_by (CAssignment) %>%
mutate (p = n / sum (n))
trt_logit <- trt_counts %>%
group_by (CAssignment) %>%
mutate (clogit = logit (cumsum (p))) %>%
group_by (out_mmrc_scale) %>%
mutate (rel_clogit = clogit - mean (clogit)) %>%
ggplot (., aes (out_mmrc_scale, rel_clogit)) +
facet_wrap ( ~ CAssignment) +
geom_point () +
labs (y = "Relative (to mean) sample cumulative logit" )
trt_logit
```
Some indication that a contrained partial proportional odds model might have better 'fit', but sample size is very small, so sticking with PO model for consistency among results.
```{r}
#| code-summary: Plot sample cumulative logits
#| fig-cap: Inspect sample cumulative logits, antiviral.
trt_counts <- fas_itt_nona_dat %>%
dplyr:: count (AAssignment, out_mmrc_scale) %>%
complete (AAssignment, out_mmrc_scale, fill = list (n = 0 )) %>%
group_by (AAssignment) %>%
mutate (p = n / sum (n))
trt_logit <- trt_counts %>%
group_by (AAssignment) %>%
mutate (clogit = logit (cumsum (p))) %>%
group_by (out_mmrc_scale) %>%
mutate (rel_clogit = clogit - mean (clogit)) %>%
ggplot (., aes (out_mmrc_scale, rel_clogit)) +
facet_wrap ( ~ AAssignment) +
geom_point () +
labs (y = "Relative (to mean) sample cumulative logit" )
trt_logit
```
# Modelling {#modelling}
## FAS-ITT
```{r}
#| label: fit-model-fas-itt
#| code-summary: Fit primary model
res <- fit_primary_model (
fas_itt_nona_dat %>% mutate (out_mmrc_scale = out_mmrc_scale + 1 ),
ordmod,
outcome = "out_mmrc_scale" ,
intercept = FALSE
)
names (res$ drws$ AOR) <- "Nafamostat"
names (res$ drws$ COR) <- intervention_labels2 ()$ CAssignment[- (1 : 2 )]
names (res$ drws$ OR) <- c ("Ineligible aspirin" , "Age \u2265 60" , "Female" , "Oxygen requirement" , "Australia/New Zealand" , "Nepal" )
```
```{r}
#| label: odds-ratio-summary-table-fas-itt
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects), FAS-ITT.
save_tex_table (
odds_ratio_summary_table (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR), "latex" ),
"outcomes/secondary/7-7-primary-model-fas-itt-summary-table" )
odds_ratio_summary_table (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR))
```
```{r}
#| label: fig-or-densities-fas-itt-model
#| code-summary: Odds ratio posterior densities
#| fig-cap: Posterior densities for odds ratio contrasts.
p <- plot_or_densities (c (res$ drws$ AOR, res$ drws$ COR)) +
labs (x = "Odds ratio (log scale)" , y = "Comparison" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" , "7-7-primary-model-fas-itt-odds-ratio-densities.pdf" )
ggsave (pth, p, width = 6 , height = 2.5 )
p
```
```{r}
#| code-summary: Odds ratio posterior summary for epoch and site
#| fig-cap: Summary of epoch and site posterior odds ratios.
p <- plot_epoch_site_terms (
res$ drws$ gamma_epoch,
res$ drws$ gamma_site,
factor (res$ dat$ region_by_site,
labels = c ("India" , "Australia \n New Zealand" , "Nepal" ))
)
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" , "7-7-primary-model-epoch-site-terms-fas-itt.pdf" )
ggsave (pth, p, width = 6 , height = 4.5 )
p
```
:::{.callout-caution collapse="true"}
### Diagnostics
```{r}
res$ fit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
### Trace plots
```{r}
mcmc_trace (res$ drws["beta" ])
mcmc_trace (res$ drws["gamma_site" ])
mcmc_trace (res$ drws["gamma_epoch" ])
```
:::
### Posterior Predictive
```{r}
#| code-summary: Posterior predictive figure
#| label: fig-fas-itt-ppc
#| fig-cap: Posterior predictive check, FAS-ITT.
y_raw <- as.integer (levels (res$ dat$ y_raw))
y_ppc <- res$ drws$ y_ppc
y_ppc_raw <- rfun (\(x) y_raw[x])(y_ppc)
ppc_dat <- bind_cols (fas_itt_nona_dat, tibble (y_ppc = y_ppc_raw))
grp_ppc1 <- function (grp) {
ppc_dat %>%
group_by (grp = !! grp) %>%
summarise (
y_1 = mean (out_mmrc_scale == 0 ),
ypp_1 = rvar_mean (y_ppc == 0 ),
y_1 = mean (out_mmrc_scale == 1 ),
ypp_1 = rvar_mean (y_ppc == 1 ),
y_2 = mean (out_mmrc_scale <= 2 ),
ypp_2 = rvar_mean (y_ppc <= 2 ),
y_3 = mean (out_mmrc_scale <= 3 ),
ypp_3 = rvar_mean (y_ppc <= 3 ),
y_4 = mean (out_mmrc_scale <= 4 ),
ypp_4 = rvar_mean (y_ppc <= 4 ),
y_5 = mean (out_mmrc_scale <= 5 ),
ypp_5 = rvar_mean (y_ppc <= 5 )
) %>%
pivot_longer (y_1: ypp_5, names_to = c ("response" , "mrc" ), names_sep = "_" , values_to = "posterior" ) %>%
mutate (mrc = as.numeric (mrc))
}
grp_ppc2 <- function (grp) {
ppc_dat %>%
group_by (grp = !! grp) %>%
summarise (
y_1 = mean (out_mmrc_scale == 0 ),
ypp_1 = rvar_mean (y_ppc == 1 ),
y_2 = mean (out_mmrc_scale == 1 ),
ypp_2 = rvar_mean (y_ppc == 2 ),
y_3 = mean (out_mmrc_scale == 2 ),
ypp_3 = rvar_mean (y_ppc == 3 ),
y_4 = mean (out_mmrc_scale == 3 ),
ypp_4 = rvar_mean (y_ppc == 4 ),
y_5 = mean (out_mmrc_scale == 4 ),
ypp_5 = rvar_mean (y_ppc == 5 ),
y_6 = mean (out_mmrc_scale == 5 ),
ypp_6 = rvar_mean (y_ppc == 6 )
) %>%
pivot_longer (y_1: ypp_6, names_to = c ("response" , "mrc" ),
names_sep = "_" , values_to = "posterior" ) %>%
mutate (mrc = as.numeric (mrc))
}
plot_grp_ppc <- function (dat) {
ggplot (dat %>% filter (response == "ypp" ), aes (x = mrc)) +
facet_wrap ( ~ grp, nrow = 1 ) +
stat_slabinterval (aes (ydist = posterior)) +
geom_point (data = dat %>% filter (response == "y" ),
aes (x = mrc, y = mean (posterior)),
colour = "red" ,
shape = 23 ) +
labs (x = "mMRC outcome" ,
y = "Probability" )
}
plot_grp_ppc <- function (dat, lab = "" , xlab = "Probability" ) {
ggplot (dat %>% filter (response == "ypp" ), aes (y = mrc)) +
facet_wrap ( ~ grp, nrow = 1 ) +
stat_interval (aes (xdist = posterior), size = 2 ) +
geom_point (data = dat %>% filter (response == "y" ),
aes (y = mrc, x = mean (posterior)),
colour = "red" ,
shape = 23 ) +
scale_x_continuous (xlab, breaks = c (0 , 0.5 ),
sec.axis = sec_axis (~ . , name = lab, breaks = NULL , labels = NULL )) +
scale_y_continuous ("WHO \n outcome" , breaks = c (1 ,3 ,5 ,7 )) +
theme (strip.text = element_text (size = rel (0.7 )),
axis.title.x = element_text (size = rel (0.7 )),
axis.text.x = element_text (size = rel (0.65 )),
axis.title.y = element_text (size = rel (0.75 )),
axis.title.x.bottom = element_blank ())
}
pp_epoch <- grp_ppc2 (sym ("epoch" )) %>%
mutate (grp = fct_inorder (factor (grp)))
pp_A <- grp_ppc2 (sym ("AAssignment" ))
pp_C <- grp_ppc2 (sym ("CAssignment" ))
pp_ctry <- grp_ppc2 (sym ("Country" ))
pp_site <- grp_ppc2 (sym ("site" )) %>%
left_join (ppc_dat %>% dplyr:: count (site, Country), by = c ("grp" = "site" ))
p0 <- plot_grp_ppc (pp_A, "Antiviral" , "" )
p1 <- plot_grp_ppc (pp_C, "Anticoagulation" , "" )
p2 <- plot_grp_ppc (pp_ctry, "Country" , "" )
p3 <- plot_grp_ppc (pp_epoch, "Epoch" , "" )
p4 <- plot_grp_ppc (pp_site %>% filter (Country == "IN" ), "Sites India" , "" )
p5 <- plot_grp_ppc (pp_site %>% filter (Country == "AU" ), "Sites Australia" , "" )
p6 <- plot_grp_ppc (pp_site %>% filter (Country == "NP" ), "Sites Nepal" , "" )
p7 <- plot_grp_ppc (pp_site %>% filter (Country == "NZ" ), "Sites New Zealand" , "" )
p <- (p0 | p1 | p2) / p3 / p4 / p5 / (p6 | p7)+
plot_layout (guides = "collect" ) &
theme (legend.position = "bottom" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" ,
"7-7-primary-model-fas-itt-ppc.pdf" )
ggsave (pth, p, width = 6 , height = 5.75 )
p
```
## ACS-ITT
```{r}
#| label: fit-model-acs-itt
#| code-summary: Fit primary model
res <- fit_primary_model (
acs_itt_nona_dat %>% mutate (out_mmrc_scale = out_mmrc_scale + 1 ),
ordmod,
outcome = "out_mmrc_scale" ,
intercept = FALSE
)
names (res$ drws$ AOR) <- "Nafamostat"
names (res$ drws$ COR) <- intervention_labels2 ()$ CAssignment[- (1 : 2 )]
names (res$ drws$ OR) <- c ("Ineligible aspirin" , "Age \u2265 60" , "Female" , "Oxygen requirement" , "Australia/New Zealand" , "Nepal" )
```
```{r}
#| label: odds-ratio-summary-table-acs-itt
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects), ACS-ITT.
save_tex_table (
odds_ratio_summary_table (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR), "latex" ),
"outcomes/secondary/7-7-primary-model-acs-itt-summary-table" )
odds_ratio_summary_table (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR))
```
```{r}
#| label: fig-or-densities-acs-itt-model
#| code-summary: Odds ratio posterior densities
#| fig-cap: Posterior densities for odds ratio contrasts, ACS-ITT.
p <- plot_or_densities (c (res$ drws$ AOR, res$ drws$ COR)) +
labs (x = "Odds ratio (log scale)" , y = "Comparison" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" , "7-7-primary-model-acs-itt-odds-ratio-densities.pdf" )
ggsave (pth, p, width = 6 , height = 2.5 )
p
```
```{r}
#| label: fig-epoch-site-acs-itt-model
#| code-summary: Odds ratio posterior summary for epoch and site
#| fig-cap: Summary of epoch and site posterior odds ratios, ACS-ITT.
p <- plot_epoch_site_terms (
res$ drws$ gamma_epoch,
res$ drws$ gamma_site,
factor (res$ dat$ region_by_site,
labels = c ("India" , "Australia \n New Zealand" , "Nepal" ))
)
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" , "7-7-primary-model-epoch-site-terms-acs-itt.pdf" )
ggsave (pth, p, width = 6 , height = 4.5 )
p
```
:::{.callout-caution collapse="true"}
### Diagnostics
```{r}
res$ fit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
### Trace plots
```{r}
mcmc_trace (res$ drws["beta" ])
mcmc_trace (res$ drws["gamma_site" ])
mcmc_trace (res$ drws["gamma_epoch" ])
```
:::
## AVS-ITT
```{r}
p1 <- avs_itt_nona_dat %>%
group_by (CAssignment = factor (CAssignment, labels = str_replace (intervention_labels ()$ CAssignment, "<br>" , " \n " ))) %>%
summarise (n = n ()) %>%
ggplot (., aes (CAssignment, n)) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Anticoagulation intervention" )
p2 <- fas_itt_nona_dat %>%
dplyr:: count (
CAssignment = factor (
CAssignment,
labels = str_replace (intervention_labels ()$ CAssignment, "<br>" , " \n " )),
out_mmrc_scale = ordered (as.integer (out_mmrc_scale), levels = 0 : 5 )
) %>%
group_by (CAssignment) %>%
mutate (p = n / sum (n)) %>%
ggplot (., aes (CAssignment, p)) +
geom_col (aes (fill = out_mmrc_scale)) +
scale_fill_viridis_d (option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = F)) +
labs (fill = "mMRC scale" , y = "Cumulative proportion" , x = "Anticoagulation intervention" )
p <- p1 | p2
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "7-7-anticoagulation-avs-itt.pdf" ), p2, height = 3 , width = 6 )
p2
```
```{r}
res <- fit_primary_model (
avs_itt_nona_dat %>% mutate (out_mmrc_scale = out_mmrc_scale + 1 ),
ordmod0,
outcome = "out_mmrc_scale" ,
vars = c ("agegte60" , "sexF" , "supp_oxy2" , "crp_tertile" ),
beta_sd_var = c (2.5 , 2.5 , 2.5 , 2.5 , 2.5 , 2.5 ),
intercept = FALSE
)
names (res$ drws$ AOR) <- "Nafamostat"
names (res$ drws$ COR) <- intervention_labels2 ()$ CAssignment[- (1 : 2 )]
names (res$ drws$ OR) <- c ("Age \u2265 60" , "Female" , "Oxygen requirement" , "CRP (2nd tertile)" , "CRP (3rd tertile)" , "CRP (unknown)" )
```
```{r}
#| label: odds-ratio-summary-table-avs-itt
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects), AVS-ITT.
save_tex_table (
odds_ratio_summary_table (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR), "latex" ),
"outcomes/secondary/7-7-primary-model-avs-itt-summary-table" )
odds_ratio_summary_table (c (res$ drws$ AOR, res$ drws$ COR, res$ drws$ OR))
```
```{r}
#| label: fig-or-densities-avs-itt-model
#| code-summary: Odds ratio posterior densities
#| fig-cap: Posterior densities for odds ratio contrasts, AVS-ITT.
p <- plot_or_densities (c (res$ drws$ AOR)) +
labs (x = "Odds ratio (log scale)" , y = "Comparison" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" , "7-7-primary-model-avs-itt-odds-ratio-densities.pdf" )
ggsave (pth, p, width = 6 , height = 2.5 )
p
```
:::{.callout-caution collapse="true"}
### Trace plots
```{r}
mcmc_trace (res$ drws["beta" ])
```
:::
# End of script.